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Abstract 

Reed instruments are modeled as self-sustained oscillators driven by the pressure inside the 
mouth of the musician. A set of nonlinear equations connects the control parameters (mouth 
pressure, lip force) to the system output, hereby considered as the mouthpiece pressure. Clarinets 
can then be studied as dynamical systems, their steady behavior being dictated uniquely by the 
values of the control parameters. Considering the resonator as a lossless straight cylinder is a 
dramatic yet common simplification that allows for simulations using nonlinear iterative maps. 

In this paper, we investigate analytically the effect of a time-varying blowing pressure on the 
behavior of this simplified clarinet model. When the control parameter varies, results from the 
so-called dynamic bifurcation theory are required to properly analyze the system. This study 
highlights the phenomenon of bifurcation delay and defines a new quantity, the dynamic oscillation 
threshold. A theoretical estimation of the dynamic oscillation threshold is proposed and compared 
with numerical simulations. 

Keywords: Musical acoustics, Clarinet-like instruments, Iterated maps, Dynamic Bifurcation, 
Bifurcation delay, Transient processes. 

1 Introduction 

One of the interests of mathematical models of musical instruments is to be able to predict certain 
characteristics of the produced sound given the gesture performed by the musician. In the case of a 
clarinet for instance, the amplitude, frequency or spectral content (the sound parameters) can be to 
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a certain extent, determined as a function of the blowing pressure and lip force applied to the reed 
(the control parameters). A basic model, such as the one introduced by Mcintyre et al. [1], allows to 
compute the amplitude of the oscillating resonator pressure from the knowledge of these two control 
parameters, giving results that follow the major tendencies observed in experiments. Several degrees of 
refinement can be added to this model, usually aiming at realistic sound and mechanical behavior. Well 
known simplifications of this model allow to study analytically the behavior of the clarinet. Simplified 
models, of course, are unable to describe or predict with refinement the exact harmonic content of 
the sound, or the influences of such important details as the reed geometry and composition or the 
vocal tract of the player. However, they can provide an understanding of the factors essential for the 
production of sound. 

The highest degree of simplification of the model (introduced in Section 2) considers a straight, 
lossless (or losses independent of frequency) resonator and the reed as an ideal spring [2, 3, 4]. With 
these assumptions, the system can be simply described by an iterated map [1]. Iterated maps often 
describe a succession of different regimes with variable periodicity. By analyzing the asymptotic values 
of these regimes it is possible to estimate: thresholds of oscillation, extinction, beating regimes, etc. 
[5], amplitudes and stability of the steady state regime [6] and phenomena of period doubling [7, 8]. 

These characteristics arise from the so-called static bifurcation theory assuming that control pa- 
rameters are constant. For example, these studies allow to find a static oscillation threshold 7^ [5] 
such that a constant regime is stable if the blowing pressure is below 7^ and a periodic regime is stable 
if it is above j s t- This behavior is static, obtained by choosing a constant blowing pressure, letting 
the system reach its final state, and repeating the procedure for other constant blowing pressures. 
Therefore, most studies using iterated map approach are restricted to a steady state analysis of the 
oscillation. They focus on the asymptotic amplitude regardless of the history of the system. 

During a note attack transient the musician varies the pressure in her/his mouth before reaching 
a quasi-constant value. During this transient the blowing pressure cannot be regarded as constant. 
In a mathematical point of view increasing the control parameter (here the blowing pressure) makes 
the system non-autonomous and results from static bifurcation theory are not sufficient to describe 
its evolution. Indeed, it is known that, when the control parameter varies, the bifurcation point - 
i.e. the value of the blowing pressure where the system begins to oscillate - can be considerably 
delayed [9, 10, 11]. Indeed, the bifurcation point is shifted from 7 s j to a larger value 7^ called dynamic 
oscillation threshold. This phenomenon called bifurcation delay is not predicted by the static theory. 
Therefore, when the control parameter varies, results from the so-called dynamic bifurcation theory 
are required to properly analyze the system. 

The purpose of this paper is to use results from dynamic bifurcation theory to describe analytically 
a simplified clarinet model taking into account a blowing pressure that varies linearly with time. In 
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particular we propose a theoretical estimation of the dynamic oscillation threshold. 

Section 2 introduces the simplified mathematical model of a clarinet and the iterated map method 
used to estimate the existence of the oscillations inside the bore of the clarinet. Some results related 
to the steady state are presented in this section. Section 3 is devoted to the study of the dynamic 
system that takes into account a time-varying blowing pressure. The phenomenon of bifurcation delay 
is demonstrated using numerical simulations. A theoretical estimation of the dynamic oscillation 
threshold is also presented and compared with numerical simulations. In Section 4 the limits of this 
approach are discussed. It is shown, when the model is simulated, that the precision (the number of 
decimal digits used by the computer) has a dramatic influence on the bifurcation delay. The influence 
of the speed at which the blowing pressure is swept is also discussed. 

2 State of the art 
2.1 Elementary model 

The model of the clarinet system used in this article follows an extreme simplification of the instrument, 
which can be found in other theoretical works [2, 4]. Although it is not suitable to describe the detailed 
content of the sound produced by the instrument, it is an useful tool to predict the magnitude of certain 
key features of the instrument, such as the threshold of oscillation, or to estimate the amplitudes of 
the pressure oscillation inside the instrument. 

This basic model separates the instrument into two functional elements. One of these is the bore, 
or resonator, a linear element where the pressure waves propagate without losses. The other is the 
reed-mouthpiece system, which is considered as a valve controlled by the pressure difference between 
the mouth and the mouthpiece. It is often called the generator and is the only nonlinear part of the 
instrument. 



Figure 1: Schematic diagram of a single-reed mouthpiece. Presentation of variables, control parameters 
and choice of axis orientation. U is the flow created by the pressure imbalance P m — P between the 
mouth and the bore, U r is the flow created by the motion of the reed, Ui n is the flow entering the 
instrument, y represents the position of the tip of the reed and H is the opening of the reed channel 
at rest. 



Reed channel 




Lip 



Mouthpiece 



3 



2.1.1 The reed-mouthpiece system 

The reed- mouthpiece system is depicted in Fig. 1. The reed is assumed to behave as an ideal spring 
characterized by its static stiffness per unit area K s . So, its response y to the pressure difference 
AP = P m — P is linear and is given by: 

From (1) we can define the static closing pressure Pm which corresponds to the lowest pressure 
that completely closes the reed channel (y = —H): 

P M = K S H. (2) 

The reed model also considers that the flow created by the motion of the reed U r is equal to zero, 
so that the only flow entering the instrument is created by the pressure imbalance between the mouth 
and the bore: 

U in = U. (3) 

The non-linearity of the reed-mouthpiece system is introduced by the Bernoulli equation which 
relates the flow U to the acoustic pressure P [12, 13]. This relation is the nonlinear characteristics of 
the exciter, given by: 



U 



AP\ /I API 

U A [1- —J J sgn(AP) if AP < P M (4a) 



L if AP > P M - (4b) 
The flow Ua is calculated using the Bernoulli theorem: 



Ua = SJ 2 ^, (5) 

where S is the opening cross section of the reed channel at rest and p the density of the air. 
Introducing the dimensionless variables and control parameters [4]: 

&p=-5- ; p = -5- > u = z c~5- ; 7 = -5- ; C = z c —, (6) 

Pm Pm Pm Pm Pm 

equation (4) becomes: 
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f C(l - 1 + P) V\l-P\ s e n (7 -P) if 7 ~P < 1 (7a) 
L if 7 -p > 1. (7b) 

The parameters 7 and (" are the control parameters of the system. An example of the function F 
is shown in Fig. 2(a). 

2.1.2 The resonator 

Assuming that only plane waves exist in the resonator and propagate linearly, the resonator can be 
characterized by its reflection function r{t). The general expression relating p(t) to u(t) through r(t) 
is: 

p(t)-u(t) = [r*(p + u)](t). (8) 

The resonator is modeled as a straight cylinder. Reflections at the open end of the resonator are 
considered perfect (no radiation losses) and viscous and thermal losses are neglected. In this case the 
reflection function becomes a simple delay with sign inversion: 



r(t) = - 

where 5 is the Dirac generalized function and r = 
with velocity c along the resonator of length I. 
With the reflection function (9), equation (8) 



-S(t-r), (9) 
: 21 /c is the round trip travel time of the sound wave 

becomes: 



p(t)-u(t) = -[ P (t-T) + U (t-T)]. (10) 

Using a discrete time formulation (the discretization is done at regular intervals r) and noting 
p{nr) = p n and u(jit) = u n , we obtain the following difference equation: 

Pn ~ u n = - (jpn-i + u n -l) ■ (11) 
2.2 Iterated map: outgoing and incoming wave representation 

In linear acoustics any planar wave can be expanded into an outgoing wave p + and an incoming wave 
p~ . Using the dimensionless variables defined in equation (6), the acoustic pressure p and flow u are 
given by: 
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P = p + + p ; u = p + — p , (12) 

Replacing in equation (11), 

p + = ^(p + u) ; p~ = -(p-u). (13) 

By combining equations (7) and (12) a nonlinear relation G between p + and p~ can be obtained: 

p+ = G{-p-). (14) 

An explicit expression of the function G was determined, for £ < 1, by Taillard et al. [8]. Fig. 2(b) 
shows an example of the function G. Using equations (12), the relation (11) becomes: 

Pn = -Pn-l- ( 15 ) 
Finally, equations (14) and (15) define the iterated map: 

p+ = G[-p-)=G{p+_ 1 ). (16) 

In the following, the variable p + will be used preferentially. The variable p can easily be calculated 
using equations (15) and (12). 




-0.5 0.5 -0.4 -0.2 0.2 0.4 

(a) function F (b) function G 

Figure 2: Nonlinear characteristics in u = F(p) representation (a) and p + = G(—p~) representation 
(b) for 7 = 0.42 and ( = 0.6. 



2.3 Results from static bifurcation theory 

The difference equation (16) can be analyzed using the static bifurcation theory, which assumes that 
the control parameters are constant. This will be hereafter referred to as the static case. The parameter 
7 will be specifically introduced as a subscript in the definition of the nonlinear characteristics (16), 
stressing that this will be the parameter of interest in the current study (£ will always consider to be 
constant): 

Vt = G, (p+_i) ■ (17) 

Some of the predictions of the static bifurcation theory that are important to this work are recalled 
in the following sections while applying them to the map of equation (17) [14, 4, 8]. 

2.3.1 Expression of the static regime and static oscillation threshold 

For all values of the control parameter 7 below a particular value of the parameter 7 called static 
oscillation threshold and noted 7^ the series p+ converges to a single value (the static regime), also 
referred to as the fixed point of G 7 . It can be found by solving the following equation: 

p+* = G 7 (p+*) . (18) 

After solving the equation we obtain : 

P + *h) = |(1~7)V7- (19) 

When the static regime is reached p+ = p^-i = ~Pn- Therefore, for the the variable p = p + 
the static regime is equal to zero. 

The static regime exists for all values of the parameter 7 but it is stable when 7 < 7^ and unstable 
when 7 > 7^. The condition of stability of the static regime [14] allowing to obtain the value of the 
static oscillation threshold is: 

|G;(p + *)|<1, (20) 
where G' is the first derivative of the function G. The value of the static oscillation threshold is finally: 

1st = \- (21) 

Beyond the oscillation threshold, other bifurcations occur, the 2-valued oscillating regime becoming 
unstable and giving rise to a 4- valued oscillating state. This cascade is the classical scenario of successive 
period doublings, leading eventually to chaos [15, 8]. The values of the parameter 7 for which appear 
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the different 2n-valued oscillating regimes depend on the value of the parameter the smaller is £, the 
earlier the 2n-valued oscillating regimes appear. When 7 = 1/2, whatever the value of £, a 2-valued 
oscillating regime reappears, the beating-reed regime . This is a particularity of model of the clarinet, 
it is due to the fact that when 7 — p > 1 (equation (7b)) the reed presses against the mouthpiece lay. 
It can be shown [4] that in this permanent regime p = ±7 (c.f. Fig. 3). 

2.3.2 Static bifurcation diagrams 

Common representations of the static bifurcation diagram for clarinets usually show the steady state 
of the pressure inside the mouthpiece p or that of its amplitude (corresponding in the lossless model 
to the absolute value of p) with respect to the control parameter 7 [5]. In this paper, calculations are 
based on p + , so that most bifurcation diagrams will represent the steady state of the outgoing wave [8]. 

Fig. 3 shows an example of these three representations of the static bifurcation diagram for £ = 0.5. 
Fig. 3 represents only the two first branches of the diagrams. The first branch corresponds to the fixed 
points of the function G 7 and the second branch represents the fixed points of the function {G~ o G~). 
On the right figure, the dashed line represents the curve of the static regime of p + , noted p + *{^). 
Oscillating regimes with higher periodicities which may appear between 7 = 1/3 and 7=1/2 are not 
represented. 




(a) (b) (c) 

Figure 3: Graphical representation of the static bifurcation diagrams for £ = 0.5. From left to right: 
diagram based on variables \p\, p and p + . 

3 Time-varying blowing pressure 

3.1 Problem statement 
3.1.1 Definitions 

Before presenting the problem, some definitions are introduced in order to avoid ambiguity in the 
vocabulary used hereafter. In the remainder of this paper, all simulations and calculations will be 
performed considering that the parameter £ is a constant, although its value may vary from one simu- 
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lation to another. The definitions presented below, used commonly in works dealing with bifurcation 
theory, can present some conflicts with that of musical acoustics. The terms that will be used in the 
remaining discussions are clarified in the following paragraphs: 



Static case The control parameter 7 is constant and the system is described by: 



Vn = G~ ( (p+-j) • 



(22) 



The steady state of the series p+ depends on the value of the control parameter 7. If 7 is smaller 
than j s t, the series tends to a static regime. To avoid confusion, the static regime will now be called 
non- oscillating static regime. If 7 is larger than 7 s t the steady state of p\ is an oscillating regime. This 
regime is called oscillating static regime. This behavior is still static, obtained by choosing a value of 
7, letting the system reach its steady state, and repeating the procedure for each value of 7. Note that, 
even if the system tends to a steady state, the initial condition p^ often induces a transient regime. 

Dynamic case The control parameter 7 is variable, now written as j n . When 7 is a linear function 
of time, the system is described by the following difference equations: 



A slowly varying parameter implies that e is arbitrarily small (e <C 1). An example of a numerical 
simulation performed on the system (23) is shown in Fig. 4 for £ = 0.5, e = 10 -3 and an initial 
condition 70 = 0.2. The initial value of the outgoing wave is p^ = G(~Pq = 0,70). Indeed, for n = 
the incoming wave p~ is clearly zero, otherwise sound would have traveled back and forth with an 
infinite velocity. 

The series p^ first shows a short oscillating transient, which will be called transient oscillating 
dynamic regime. This oscillation decays into a non-oscillating dynamic regime. Beyond a certain 
threshold, a new oscillation grows, giving rise to the final oscillating dynamic regime. 




(23b) 
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n n 
(a) (b) 

Figure 4: Numerical simulation performed on the system (23). (a) complete orbit of the series and (b) 
zoom near the non-oscillation dynamic regime. Q = 0.5, e = 10~ 3 , 70 = 0.2 and p^ = G(0, 70). 

This paper will focus on the transition (i.e. the bifurcation) from the non-oscillating dynamic 
regime to the final oscillating dynamic regime. The value of the parameter 7 for which the bifurcation 
occurs is called dynamic oscillation threshold, noted 7^. 

3.1.2 Bifurcation delay 

Bifurcation delay occurs in nonlinear-systems with time varying control parameters. Fruchard and 
Schafke [11] published an overview of the problem of bifurcation delay. 

In fig. 5, the system (23) was simulated numerically, showing the time evolution of the series p^ and 
of the control parameter 7„ (cf. Fig. 5(a)). To better understand the consequence of a time- varying 
parameter, the orbit of the series p^ is plotted as a function of the parameter 7 n - in this case the 
evolution of the system can be interpreted as a dynamic bifurcation diagram. This is compared to the 
static bifurcation diagram in Fig. 5(b). We can observe that the static and the dynamic oscillation 
diagrams coincide far from the static oscillation threshold 7^. However, in the dynamic case, we can 
see that the orbit continues to follow closely the branch of the fixed point of function G throughout a 
remarkable extent of its unstable range, i.e. after 7^: the bifurcation point is shifted from the static 
oscillation threshold 7^ to the dynamic oscillation threshold 7^. The term bifurcation delay is used to 
state the fact that the static oscillation threshold j s t is smaller than the dynamic oscillation threshold 
Jdt- 
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(a) (b) 



Figure 5: (a) Time evolution of the series p+ and of the control parameter 7 n .(b) Comparison between 
the series p^ and the static bifurcation diagram as a function of 7 n . £ = 0.5, e = 10~ 4 , 70 = and 
p+ = G(0, 7 o). 

Non-standard analysis has been used in the past to study the phenomenon of bifurcation delay 
[16, 17], explaing that one of the causes of the bifurcation delay is the exponential proximity between 
the orbit of the series p+ and the curve of the the fixed point of G. Other studies of bifurcation delay 
using standard mathematical tools - mathematics [10, 18] or physics publications [9, 19] - explain 
bifurcation delay as an accumulation of stability during the range of 7 for which the fixed point of G 
is stable (i.e. < 7„ < 7^). The dynamic oscillation threshold therefore appears as the value of the 
parameter 7 at which the stability previously accumulated is compensated. 

In musical acoustics literature some papers present results showing the phenomenon of bifurcation 
delay without never making a connection to the concept of dynamic bifurcation. For example this 
phenomenon is observed in simulations of clarinet-like systems using a slightly more sophisticated 
clarinet model (Raman's model) [20]. Raman's model takes losses into account although they are 
assumed to be independent of frequency (see [5] for further explanation). Bifurcation delay can also 
explain the difficulty in estimating the static oscillation threshold by using a slowly variable blowing 
pressure [21]. In a preliminary work [22], bifurcation delays were experimentally observed in a clarinet- 
like instrument. 

3.2 Analytical study of the dynamic case 

This section presents an analytical description of a clarinet-like system in a dynamic case. The notion 
of invariant curve (0(7,e)), invariant under the mapping (23), will be needed for this study. The 
study of the stability of the invariant curve allows to define an analytical estimation of the dynamic 
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oscillation threshold. A generic method to calculate the invariant curve is given by Baesens [10] a , based 
on a perturbation method [23]. 

3.2.1 Invariant curve 

The invariant curve ^(7, e) is invariant under the mapping (23), satisfying the following equation: 

0(7,e) = G(0( 7 -e,e), 7 ). (24) 

This curve plays a similar role for the dynamic system as fixed points for the static system, attracting 
or repelling the orbits. It is independent of the initial condition. 

First of all, the invariant curve is expanded into a power series of e, here truncated to the first 
order: 



0(7, e) « (£0(7) +^1(7). 



(25) 



Fig. 5 shows that, during the dynamic phase, the orbit of the series p£ closely follows the curve of 
the fixed points of G. This allows to linearize function G around the curve of the fixed points p + *('j): 



using the notation 



G(x, 7 ) ~ G(p + *( 7 ),7) + [x-p + *(l)] d x G(p+*(j),j) , 

dG{x,y) 



d x G(x,y) 



dx 



(26) 



(27) 



and knowing that G(p + *('y),j) = p + *{"f) (cf. equation (18)). Finally, using a Taylor expansion of 
0(7 — e, e) equation (24) is successively solved for the functions ^0(7) and ^1(7), yielding: 



<K7, e) ~P + *(l) + e 



dp+*(j) d x G{p+*{ 1 ), 1 ) 



d 7 d x G(p+*( 1 ), 1 )-l 
Using the explicit expressions of p + * and dp + * /dj we have: 



(28) 



<M7,e) 



(1-7)^7 



(3 7 - 1) d x G(p+*( 7 ), 7 ) 



(29) 



2^7 d x G(p+*( 1 ), 1 )-l_ 
More details about the calculation of the invariant curve are given in B. 

To simplify the notation, in the rest of the document the invariant curve will be noted ^(7). Its 
dependency on parameter e is not explicitly stated. 

a In [10] the invariant curve is called adiabatic invariant manifold. 



12 



3.2.2 Stability of the invariant curve and theoretical estimation of the dynamic oscilla- 
tion threshold 

A theoretical estimation of the dynamic oscillation threshold is done by identifying the value of 7 for 
which the invariant curve looses its stability. The invariant curve is said to be unstable when the orbit 
of of the series escapes from the neighborhood of the invariant curve ^(7, e). 

To investigate the stability of the invariant curve ^(7, e), the function G in equation (23a) is 
expanded in a first-order Taylor series around the invariant curve [10]: 

Pn = G{p+_ x ,~f n ) 

« G (0( 7n - e), 7 n) + [pt-i ~ <M7n - e)] d x G (0( 7 „ - e), ln ) . (30) 

A new variable is defined that describes the distance between the actual orbit and the invariant 
curve: 



w n =pt ~ Hln), (31) 
and using equation (24), equation (30) becomes: 



w n = w n -id x G (4>{^ n - e),7 n ) . (32) 
The solution of equation (32) is formally: 

n 

w n = w \\d x G((j)(ji - e) )7i ) , (33) 

8=1 

for n > 1 and where wo is the initial value of w n . The absolute value of w n can be written as follow: 

\w n \ = \wo\ exp ^Y]]n\d x G (0( 7i - e), 7i )|^ . (34) 
Finally, using Euler's approximation the sum is replaced by an integral: 

\w n \ ~ \w \ exp ( I In \d x G (0(7 - e), V) | — V (35) 

Equation (35) shows that the variable p + starts to diverge from the invariant curve 0(7, e) when 
the argument of the exponential function changes from negative to positive. Therefore, the analytical 
estimation of the dynamic oscillation threshold 7^ is defined by: 
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/ ' hi\d x G(<j>( 1 , -e), 1 , )\dh , = O t (36) 

J 70+e 

where 70 is the initial value of 7. This result can be deduced from [10] (equation (2.18)), it may also 
be obtained in the framework of non-standard analysis [24]. 

The theoretical estimation 7^ of the dynamic oscillation threshold depends on the initial condition 
70 and on the increase rate e, it is therefore written 7^(70, e)- 

A numerical solution 7^(70, e) of the implicit equation (36) is plotted in Fig. 6 as a function of the 
initial condition 70 and for e = 10 -2 , 10 -3 and 10 -4 . 7^+ can be much larger than static oscillation 
threshold 7^ = 1/3 for small initial conditions 70. When the initial condition value 70 increases, 7^ 
approaches the static threshold. Fig. 6(d) shows that the bifurcation delay seems to be independent 
of the increase rate e if this value is sufficiently small (typically < 10 -3 ). 

Equation (36) states that when 7 = 7^ we have \w n \ ~ \w$\, providing a good estimation of the 
dynamic oscillation threshold 7^ if |^o| is sufficiently small, i.e. if p£ is sufficiently close to (^(70). 
70 = can be problematic since 0(0, e) = —00, but a single iteration is sufficient to bring the orbit to 
a neighborhood of the invariant curve. Therefore, we make the assumption that 

T$(0,e)«T$(e,e). (37) 

A non-exhaustive study done by running a few simulations shows that for e = 10 -4 the error in 7^ 
due to this approximation is under 10 -8 , rising to 10 -7 when e = 10 -3 and 2 x 10 -5 when e = 10~ 2 . 

3.3 Benchmark of theoretical estimators for the dynamic threshold 

Multiple criteria can be associated to the beginning of the oscillating regime. For instance, the os- 
cillations can start before the series departs from the vicinity of the invariant curve as described in 
equation (36). Moreover, because of the approximation used between equations (35) and (36), the 
value of 7 = 7^ may not be an accurate estimation of the value at which the orbit departs from this 
vicinity. 

For comparison, a dynamic oscillation threshold (noted J^t™) is calculated by simulating system 
(23) and compared with 7$^. When the orbit of the series p+ is periodic, the sign of the second order 
difference of changes sign at each iteration (i.e. the curve of changes from upward to downward 
concave). In discrete time formulation the second order difference is given by: 

(dV), = (pf - Pti) - (Pt-i - PU) ■ ( 38 ) 
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Therefore, 7^ um , the oscillation threshold measured in numerical simulations, is reached when 

(dV)^ (d 2 P + ) l < 0, (39) 

is satisfied. 




5-1CT 2 0.1 0.15 0.2 0.25 0.3 5 . i(r 2 0.1 0.15 0.2 0.25 0.3 

7o 7o 
(c) e = 10- 4 (d) 



Figure 6: Plot of j^t as a function of the initial condition 70, for different values of the slope e. (a), (b) 
and (c): solid lines are the theoretical prediction rfh calculated from equation (36). Gray " *" markers 
represent the value 7^" m for which the system begins to oscillate, (d): combination of the previous 
theoretical predictions. represent the highest 70 for which the system has enough time to reach a 
non-oscillating dynamic regime. 

Then, in Fig. 6, 7^" m is compared with 7^ (gray " * " markers). In some cases the series never 
reaches the non-oscillating dynamic regime. An example of such situations is shown in Fig. 7. The 
values of 7^" m corresponding to the last initial values 70 for which the system has enough time to 
reach the non-oscillating dynamic regime are circled. 

Fig. 6 shows that for e = 10 -4 the theoretical result 7^ provides a good estimation of the observed 
dynamic oscillation threshold. For e = 10 -3 , the theoretical estimation is also good if the the initial 
condition is sufficiently small but as 70 gets closer to the static threshold 7^ the system begins to 
oscillate before 7 = 7S+. Finally, for e = 10 -2 , 7J"" 1 is always smaller than 7^. 
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Figure 7: Representation of the series as a function of 7„ for Q = 0.5, e = 10~ 3 , 70 = 0.3 and 
p+ = G(0, 7o ). 

4 Limit of the model: influence of the precision 

The phenomenon of bifurcation delay is very sensitive to noise: either numerical noise (round-off errors 
of the computer) or experimental noise (due to turbulence for instance). Indeed, when the static 
threshold is exceeded the system is very unstable. As a result, to observe bifurcation delay with 
numerical simulations and compare to theoretical results, it is necessary to perform calculations using 
a very high precision, as was done previously in this paper. For lower precisions the bifurcation delay 
can be considerably reduced (see [17] for an example in the logistic map). 

Fig. 8 shows an example of numerical simulation performed on system (23). Fig. 8(a) and Fig. 8(b) 
differ only in the numerical precision (i.e. the number of decimal digits) used to calculate the orbit. 
The choice of the precision is possible using mpmath, the arbitrary precision library of Python. Fig. 8(a) 
was obtained using a precision of 5000 decimal digits, in this case 7^? gives a good estimation of the 
bifurcation point, as it has already been shown in Fig. 6. On the other hand, using a precision of 15 
decimal digits (Fig. 8(b)), the bifurcation delay is considerably reduced and the theoretical estimation 
of the dynamic oscillation threshold is not valid. 
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Figure 8: Representation of the series for Q = 0.5, e = 10 4 , 70 = 0, p$ = G(0,7o) and for two 
different values of the precision. 

To highlight the influence of the precision 7™ m is calculated for different precisions. Results are 
plotted in Fig. 9 and compared to the analytical values 7^ and 7^. 
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Figure 9: Graphical representation of 7™ m for different precisions (prec. = 7, 15, 100, 500 and 5000) 
and for e = 10~ 4 . Results are also compared to analytical static and dynamic thresholds: j s t and 7^. 
C = 0.5 and 70 = 0. 



The first thing to observe in Fig. 9 is the very high sensitivity of 7^" m to precision, yet all the values 
of 7^" m appear between 7^ and 7^. For the lowest precision (7 decimal digits) the bifurcation delay 
disappears and 7^" m = 7^. If the precision is very high (5000 decimals) 7^" m = 7^. Therefore, 7^ 
can be interpreted as the limit of the bifurcation delay when precision tends to infinity. In cases with 
intermediate precisions (prec. = 15, 100 and 500) the bifurcation delay increases with the precision. 

The sensitivity to the precision depends on the value of the increase rate e: Fig. 10 plots 7j" m with 
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respect to e for different values of the numerical precision. Results are also compared with 7^ and 7^. 
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e 

Figure 10: Graphical representation of 7^" m as a function of e for £ = 0.5 and using five different 
precisions. A logarithmic scale is used in abscissa. 

As above, for the lowest precision (7 decimals) the bifurcation delay disappears when e is sufficiently 
small. Indeed, r f^ m is constant and equal to 7^. Then bifurcation delay occurs and increases with 
e. The case of the highest precision (5000 decimals) is identical to an analytical case which would 
correspond to infinite precision. When e is sufficiently small, the curves of 7^" m an d 7^ overlap. 
In this case when e is small 7^" m is almost constant suggesting that the bifurcation delay does not 
depend on the increase rate, as previously shown in Fig. 6(d). Then, still in the case of a precision 
of 5000 decimals, 7J"" 1 decreases for increasing e, and 7*3? also decreases but to a lesser extent. For 
intermediate precisions (15, 100 and 500 decimals) the curve of 7^ f ?tm first increases before stabilizing 
close to the curve of 7^. For a given value of the precision, the larger the e, the smaller is the 
accumulation of round-off errors created by the computer to reach a certain value of 7. This explains 
why the bifurcation delay first increases if the precision is not sufficiently high to simulate an analytic 
case. Beyond a certain value, all curves coincide with the one corresponding to the highest precision. 
That means that the system has reached the pair of parameters [precison ; e] needed to simulate an 
analytic case. 

5 Conclusion 

When considering mathematical models of musical instruments, oscillation threshold obtained through 
a static bifurcation analysis may be possibly very different from the threshold detected on a numerical 
simulation of this model. 

For the first time for musical instruments, the differences between these two thresholds have been 
interpreted as the appearance of the phenomenon of bifurcation delay in connection with the concept 
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of dynamic bifurcation. 

Theoretical estimations of the dynamic bifurcation provided in this paper have to be compared 
with care to numerical simulations since the numerical precision used in computations plays a key role: 
for numerical precisions close to standard machine precision, the bifurcation towards the oscillating 
regime can occur at significantly lower mouth pressure values (while different most of the time from 
the threshold obtained through static bifurcation theory). Moreover, in that case, the threshold at 
which the oscillations start becomes more dependent on the increase rate of the mouth pressure. 

The dependency on precision can be linked to the influence of noise generated by turbulence as the 
musician blows into the instrument. This would explain why the delays observed in artificially blown 
instruments are shorter than the predicted theoretical ones. This will be the subject of further work 
on this subject, as well as the validity of these results for smoother curves of variation of the mouth 
pressure. 

Moreover, in the light of results presented here for a basic model of wind instruments, varying 
the blowing pressure (even slowly) does not appear as the best way to experimentally determine Hopf 
bifurcations (static). In a musical context, since the blowing pressure varies through time, the dynamic 
threshold is likely to give more relevant informations than the static threshold, even if, in a real situation 
the influence of noise must be considered. 

As a final remark, the simplistic model used in this work only describes one point per half-period 
of the sound played by the instrument. It is thus not suitable to describe different regimes (whose 
frequencies are harmonics of the fundamental one) that can be obtained by the instrument. However 
a simple extension of this model calculating the orbits of different instants within the half-period may 
be able to provide some insight on this subject. 
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A Table of notation 



Physical variables 



Symbol 


Explanation 


Unit 


Zc 


characteristic impedance 


Pa-s-m -1 


K s 


static stiffness of the reed 


Pa-m -1 


Pm 


static closing pressure of the reed 


Pa 


H 


opening height of the reed channel at rest 


m 


U 


flow created by the pressure imbalance between the mouth and the 


111 -s 




mouthpiece 




U r 


flow created by the motion of the reed 


s —1 
m -s 




flow at the entrance of the resonator 


s —1 
m -s 


U A 


flow amplitude parameter 


m 3. g -l 


p 


musician mouth pressure 


Pa 


P 


pressure inside the mouthpiece 


Pa 


AP 


pressure difference P m — P 


Pa 


y 


displacement of the tip of the reed 


m 


T 


round trip travel time of a wave along the resonator 


m 
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DlMENSIONLESS VARIABLES 



Symbol Associated physical variable 



7 musician mouth pressure 

C flow amplitude parameter 

u flow at the entrance of the resonator 

p pressure inside the mouthpiece 

r reflexion function of the resonator 

p + outgoing wave 

p~ incoming wave 

p + * non-oscillating static regime of p + (fixed points of the function G) 

(p invariant curve 

w difference between p + and <f> 

e increase rate of the parameter 7 

7si static oscillation threshold 

7dt dynamic oscillation threshold 

7^ theoretical estimation of the dynamic oscillation threshold 

7™ m value of 7 when the system begins to oscillate (calculated numerically) 



Nonlinear characteristic of the embouchure 



Function 


Associated representation 


Definition 


F 


{u;p} 


u = F(p) 


G 


{p + ; p-} 


p+ = G(-p-) 


Invariant 


curve 





The invariant curve 1^(7, e) is invariant under the mapping (23), it therefore satisfies the following 
equation: 

#y,e) = G(0(7-e,e),7). (40) 

First of all, the invariant curve is expanded into a power series of e and only he first-order is 
retained: 
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0(7, e) « 0o(7) + e0i(7). 
Secondly, the function G is linearized around the curve J> + *( 7 ) of the fixed points: 



(41) 



G(x,i) « G(p+*(7),7) + [^-p + *(7)]9 :!; G(p + *( 7 ), 7 ) 

= p +*( 7 ) + [ a; -p+*( 7 )]a :l .G(p + *( 7 ), 7 ) 



(42) 
(43) 



where 



d x G(x,y) 



9G(x,y) 
dx 



Then, we make a Taylor expansion of 0( 7 — e, e): 



0(7 -e,€) « 0( 7 ,e)-e^( 7 ,e) + O(6 2 ); 

= 0o(7) + ^i(7)- e ^^ + O(6 2 ) 



Finally, neglecting the second-order terms in e, equation (40) becomes: 



0o( 7 ) + e0i(7) =P + *(7) + 



0o(7) + ^i(7)-e^^-P + *(7) 



(44) 



(45) 
(46) 



a,G(p+*(7), 7 ). (47) 



To obtain the approximate analytical expression of the invariant cure 4>, equation (47) is successively 
solved for the functions ^0(7) anci 4>\ (7)- 
To order 0, we have to solve: 



0o(7) =P + * (7) + [0o(7) ~P + * (7)] d x G (p + *( 7 ), 7 ) • 
Therefore the expression of ^0(7) is: 



(48) 



<£o (7) 



P + *(7) 
P + *(7). 



l-^G(p+*(7),7) 
l-^G(p+*( 7 ),7)' 



(49) 
(50) 



To order 1, we have to solve: 
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Mi) = 



*.M- 8 * ,(7) 



01 (7) 



c?7 



d x G{p + *( 7 ), 7 ); (51) 
^G(p+*( 7 ),7), (52) 



and therefore: 



0l(7)= a 7 ^(p + * (7 ),7)-i (53) 



Finally the expression of the invariant curve is: 
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